%% Sample Input for Coarse Mesh Diffusion Solver COMEDS2D
% !! This input will correspond to a version which does
%    not use interpolation, generating new responses for
%    each outer iteration !!
% j.roberts 2/23/2010
clear;
%profile on
name = 'sample_input_v1';

%--------------------------------------------------------------------------
% BASIC PROBLEM INPUT VARIABLES

% number of element types
numtypes = 4;
% number of groups
numg  = 2; 
% number of materials
numm = 4;
% legendre order
order = 2;

%--------------------------------------------------------------------------
%  READ IN DATA -- HERE, JUST DEFINED 
buck = 0.8e-4;
%        D      sigR    nu*sig        scattering
dat = [ 1.500   0.030   0.000  0.000  0.000  
        0.400   0.080   0.135  0.020  0.000
        1.500   0.030   0.000  0.000  0.000
        0.400   0.085   0.135  0.020  0.000
        1.500   0.030   0.000  0.000  0.000  
        0.400   0.130   0.135  0.020  0.000 
        2.000   0.040   0.000  0.000  0.000
        0.300   0.010   0.000  0.040  0.000 ];     
dat(:,2)=dat(:,2)+buck*dat(:,1);    
% absorption cross-section    
abdat = [0.010
         0.080
         0.010
         0.085
         0.010
         0.130
         0.000
         0.010]+buck*dat(:,1);
     
%--------------------------------------------------------------------------     
%  PLACE THE ELEMENTS
elements = [ 3 2 2 2 2 2 2 3 3 2 2 2 2 1 1 4 4   % 1
             2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 4 4   % 2
             2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 4 4   % 3
             2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 4 4   % 4
             2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 4 4   % 5
             2 2 2 2 2 2 2 2 2 2 2 1 1 4 4 4 4   % 6
             2 2 2 2 2 2 2 2 2 2 2 1 1 4 4 4 4   % 7           
             3 2 2 2 2 2 2 3 3 1 1 1 1 4 4 0 0   % 8
             3 2 2 2 2 2 2 3 3 1 1 1 1 4 4 0 0   % 9
             2 2 2 2 2 2 2 1 1 1 1 4 4 4 4 0 0   % 10
             2 2 2 2 2 2 2 1 1 1 1 4 4 4 4 0 0   % 11
             2 2 2 1 1 1 1 1 1 4 4 4 4 0 0 0 0   % 12
             2 2 2 1 1 1 1 1 1 4 4 4 4 0 0 0 0   % 13       
             1 1 1 1 1 4 4 4 4 4 4 0 0 0 0 0 0   % 14
             1 1 1 1 1 4 4 4 4 4 4 0 0 0 0 0 0   % 15
             4 4 4 4 4 4 4 0 0 0 0 0 0 0 0 0 0   % 16
             4 4 4 4 4 4 4 0 0 0 0 0 0 0 0 0 0   % 17
             ];  
% elements =[  3  2  2  2  3  2  2  1  4  
%              2  2  2  2  2  2  2  1  4
%              2  2  2  2  2  2  1  1  4
%              2  2  2  2  2  2  1  4  4 
%              3  2  2  2  3  1  1  4  0
%              2  2  2  2  1  1  4  4  0
%              2  2  1  1  1  4  4  0  0
%              1  1  1  4  4  4  0  0  0
%              4  4  4  4  0  0  0  0  0 ];          
% 1.105835747253969

if (max(max(elements))>numtypes)
    disp('numtypes < largest value of elements')
    return
end

%  Note about numbering: 
%            1 3 5
%            2 4 6 ...etc... bottom to top, since array is sideways

%--------------------------------------------------------------------------
% SET THE BOUNDARY CONDITIONS:  0==vacuum, 1==reflective
%  bc(1) --> left,   bc(2) --> right
%  bc(3) --> bottom, bc(4) --> top
bc(1) = 1;
bc(2) = 0;
bc(3) = 1;
bc(4) = 0;

% modularizing:
%  run geomsetup
%  run solver
%  run output

%--------------------------------------------------------------------------
% SINGLE ELEMENT DISCRETIZATION
% -- this does not necessarily need to be the same for all elements
%    but for now, it's easiest to work with squares
Y=40;
xcm  = [   0     10  ];
xfm  = [     1       ]*Y;
ycm  = [   0     10  ];
yfm  = [     1       ]*Y;
% compute number in each region
numycm = length(yfm); % number of y regions
numxcm = length(xfm); % number of x regions
mt = zeros(numxcm,numycm,numtypes);

% ELEMENT 1 - material specification
mt(:,:,1)   = [  1 ]; 
% ELEMENT 2 - material specification
mt(:,:,2)   = [  2  ]; 
% ELEMENT 3 - material specification
mt(:,:,3)   = [  3 ]; 
% ELEMENT 4 - material specification
mt(:,:,4)   = [  4  ]; 
                                      
tol     = [1.e-6 1.e-6]; % tolerance for norms             
maxit   = 15;

swpi = 1;
swnt = 1;
swgm = 0;

for order=3:3
input = makeinput( name, numtypes, numg, numm, order, dat, abdat, ...
                   elements, bc, xcm, xfm, ycm, yfm, mt, tol, maxit, ... 
                   swpi, swnt, swgm ) ; 
              
% return  
%-------------------------------------------------------------------------
% RUN SOLVER and OUTPUT EITHER TO FILE OR SCREEN      

% tic
% [x1,normk1,it,nrm1] = solverPI(input);
% fprintf(1,'  final keff = %12.8f\n',x1(end-1));
% toc

tic
input.maxit=1;
[xseed,tmp,tmp,tmp]=solverPI(input);
[x2,normk2,it2,nrm2] = solverNT(input,xseed);
fprintf(1,'  final keff = %12.8f\n',x2(end-1));
toc
end
%twoDtwo_sparse_mod
%profile viewer

% possible solver runtime parameters
%     max outers (for iterative)
%     tolerances
%     precondition
%     which solver(s) to use
%

% single element, reflect on 4
%   Y    O     k1                 k2                  %      
%   6    3     1.083028163054769  1.082750633028797   0.025631943081439
%  12    3     1.079839200602866  1.079703240094721   0.012592396048852
% single element, vac on right and top
%   Y    O     k1                 k2                  %      
%   6    3     0.625578053601234  0.625658197825914   -0.012809
%  12    3     0.630559906657183  0.630619579976003   -0.009462
%  two
%  Y    o     k1                 k2                  %     
%  2    3     1.02473269         1.021663636123922  0.300397
%  4    3     1.01584335         1.015112556936627  0.071991
%  8    1     1.01349115         1.013164642941905  0.032226
%twoDtwo_sparse_mod
